#!/usr/bin/perl

# (c) Michael Hiller, 2007

use strict;
use warnings;
use Cwd;
use Data::Dumper;

# call perl ComputePUvaluesSlidingWindows.perl -s sequence [-m motiflen]
#
# program computes the local PU values for all n-mers using RNAfold where n is the motiflen (n=6 per default)
# different (symmetrical) context lengths are used and values are averaged (default length 11 .. 30 nt)
# output consists of a gif file that contains a graphical representation of the PU values (each data point represents the substring starting at this position)


my $ViennaRNA = "/usr/local/bin/"; #"/users/bioinfo/zhouyu/bin/bin/"; #"./";
my $KT = 0.616320775;


###################################################################################################
# compute prob that the given region is unpaired using RNAup
###################################################################################################
sub getUnpairedProb {
	my ($seq,$motiflen,$motifPosition) = ($_[0],$_[1],$_[2]);
	
	my @res = readpipe("$ViennaRNA/RNAup $seq -u $motiflen -d 2");

	# values for motifs are given at the end position of the motif --> get the end position
	$motifPosition += $motiflen;		# 1-based   (on a 0-based scheme the end pos is start+length-1  --> add 1 to get the pos on a 1-based schema

	my $line = $res[$motifPosition];		# perl arrays begin with 0 but we skip the first header line --> direct access with the motifPosition
	my @f = split(/\s+/, $line);
	
	return $f[1];
}



###################################################################################################
# compute prob that the given region is unpaired using RNAfold
###################################################################################################
sub getUnpairedProbRNAfold {
	my ($seq,$motiflen,$motifPosition) = ($_[0],$_[1],$_[2]);

	# fold with constraint that the motif is unpaired
	my $constraint = "." x length($seq);						# x = unpaired, . = no constraint
	for (my $i=$motifPosition; $i<$motifPosition+$motiflen; $i++) {
		substr($constraint, $i,1) = "x";
	}
	
	my @res = readpipe "echo \"$seq\n$constraint\" | $ViennaRNA/RNAfold -d 2 -p0 -C ";
#	print "$seq\n$constraint\n";
#	print "@res\n";
	
	my $energy_constraint = 100116;		# just init
	if ($res[2] !~/ free energy of ensemble = (.*[-\d\.e]+)\s+kcal\/mol.*/) {
		print "INTERNAL ERROR: cannot parse ensemble from $res[2]<br>\n";
		exit -1;
	}else{
#		print "$res[2] ---> $1\n";
		$energy_constraint = $1;	
	}

	# fold without constraints
	@res = readpipe "echo \"$seq\" | $ViennaRNA/RNAfold -d 2 -p0 ";
	my $energy_all = 100116;		# just init
	if ($res[2] !~/ free energy of ensemble = (.*[-\d\.e]+)\s+kcal\/mol.*/) {
		print "INTERNAL ERROR: cannot parse ensemble from $res[2]<br>\n";
		exit -1;
	}else{
#		print "$res[2] ---> $1\n";
		$energy_all = $1;	
	}

#	print "****** $energy_all  $energy_constraint\n";
	if (($energy_constraint == 100116) || ($energy_all == 100116)) {
		return "-";
	}
	
	# compute prob of set of all structures with unpaired constraints (= Z_unpaired / Z; Z = partition function)
	my $prob_unpaired = exp( ((-1*$energy_constraint) - (-1*$energy_all)) / $KT);
#	print "******        $prob_unpaired\n";
	return $prob_unpaired;
}




#########################################################################################################
# compute average of vector
#########################################################################################################
sub getAve {
	my @vector = @_;
	if ($#vector < 0) {
		return 0;
	}

	my $sum = 0;
	my $count = 0;	
	for (my $i=0; $i<=$#vector; $i++) {
		$sum += $vector[$i];
		$count++;
	}
	$sum /= $count;
	return $sum;
}




###################################################################################################
# 														main
###################################################################################################

my $minContextLen=11;
my $maxContextLen=30;
my $inputsequence = "";
my $motiflen = 6;
my $FilenameHeader = "";

# determine input parameters
for (my $i = 0; $i < @ARGV; $i++){
	if ($ARGV[$i] eq "-s") {
		if (defined($ARGV[$i+1])){
			$inputsequence = $ARGV[$i+1];
		}
	}
	if ($ARGV[$i] eq "-m") {
		if (defined($ARGV[$i+1])){
			$motiflen = $ARGV[$i+1];
		}
	}
	if ($ARGV[$i] eq "-c") {
		if (defined($ARGV[$i+1])){
			$minContextLen = $ARGV[$i+1];
		}
		if (defined($ARGV[$i+2])){
			$maxContextLen = $ARGV[$i+2];
		}
	}
	if ($ARGV[$i] eq "-H") {
		if (defined($ARGV[$i+1])){
			$FilenameHeader = $ARGV[$i+1];
		}
	}
}




# check for input of context parameters and filename
if (($minContextLen <= 0) || ($maxContextLen <= 0)) {
	print "Wrong context parameters given: -c ContextLengthStart ContextLengthEnd\n";
	exit(-1);
}
if ($inputsequence eq ""){
	print "You have to specify an input sequence: -s seq\n";
	exit(-1);
}
if ($FilenameHeader eq ""){
	print "INTERNAL ERROR: no filename header given<br>";
	exit(-20);
}

if ($FilenameHeader !~ /.*\/$/) {
	$FilenameHeader .= "/";
}
#print "working with $inputsequence  motiflen: $motiflen  context-lengths: $minContextLen .. $maxContextLen   Header $FilenameHeader <br>\n";



# print the sequence into a file
my $outfile = $FilenameHeader."seq.fa";
open OUTFILE, ">$outfile" || die "cannot open $outfile\n";
print OUTFILE ">sequence\n$inputsequence\n";
close OUTFILE;


# stores PU values for all motifs
my @PUvalues;

# for all motifs of length motiflen
for (my $start=$maxContextLen; $start<=length($inputsequence)-$maxContextLen-$motiflen; $start++) {
	my $up = substr($inputsequence,$start-$maxContextLen, $maxContextLen);
	my $motif = substr($inputsequence,$start,$motiflen);
	my $down = substr($inputsequence,$start+$motiflen, $maxContextLen);
#	print "current motif start: $start  \n$up $motif $down\n";

	my @curValues;

	# do it for all context lengths 
	for (my $curContextLen = $minContextLen; $curContextLen <= $maxContextLen; $curContextLen++) {

		my $curUp = substr($up,length($up)-$curContextLen,$curContextLen);
		my $curDown = substr($down,0,$curContextLen);
		my $motifPos = length($curUp);

		# check if sequences are long enough
		if ((length($curUp) < $curContextLen) || (length($curDown) < $curContextLen)) {
			print "ERROR: seq_up or seq_down smaller than current context length $curContextLen  ($curUp  $curDown)\n";
			exit(-1);
		}

		# compute prob unpaired using RNAup
		my $PU = getUnpairedProbRNAfold($curUp.$motif.$curDown,$motiflen, $motifPos);
#		my $PU = getUnpairedProb($curUp.$motif.$curDown,$motiflen, $motifPos);
		# in case RNAfold with constraints fails, compute PU with RNAup
		if ($PU eq "-") {
			$PU = getUnpairedProb($curUp.$motif.$curDown,$motiflen, $motifPos);
		}
		push @curValues, $PU if ($PU ne "-");		
#		print "PU for contextlength $curContextLen: $PU\n";
	}
	
	my $avePU = getAve(@curValues);
#	print "FINAL average PU: $avePU\n\n";
	push @PUvalues, $avePU;
}

#print Dumper(@PUvalues);

# print the data values
my $index = 0;
for (my $start=$maxContextLen; $start<=length($inputsequence)-$maxContextLen-$motiflen; $start++) {
	printf "%d %6.4f %s\n", $start+1, $PUvalues[$index], substr($inputsequence, $start, $motiflen);
	$index++;
}

